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Abstract 

We use the single-cluster Monte Carlo update algorithm to simulate 
the three-dimensional classical Heisenberg model in the critical region 
on simple cubic lattices of size with L = 12, 16, 20, 24, 32, 40, and 
48. By means of finite-size scaling analyses we compute high-precision 
estimates of the critical temperature and the critical exponents, using 
extensively histogram reweighting and optimization techniques. Mea- 
surements of the autocorrelation time show the expected reduction 
of critical slowing down at the phase transition. This allows simula- 
tions on significantly larger lattices than in previous studies and con- 
sequently a better control over systematic errors in finite-size scaling 
analyses. 
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The critical behaviour of the three-dimensional (3D) classical Heisenberg 
model, as one of the simplest spin models, has been investigated by a variety 
of approaches. Despite this fact there are still some discrepancies to be 
resolved. Motivated by conflicting estimates of its critical coupling (3c on 
a simple cubic lattice coming from widely accepted high-temperature series 
expansion analyses {(3c ~ 0.6916) and more recent transfer-matrix (TM) 
Monte Carlo (MC) investigations [| {(3c = 0.6922(1) and (3c = 0.6925(3)), 
Peczak et al. (PFL) have recently undertaken a high statistics MC study 
of this model on cubic lattices of sizes up to = = 24^. By simulating 
the system with the standard Metropolis algorithm and making extensive 
use of multi-histogram techniques they could not decide between the 
two alternatives. Rather, from a finite-size scaling (FSS) |0 analysis of the 
crossing points of Binder's cumulant they claimed an even larger value of 
(3c = 0.6929(1). Moreover, extrapolations of the locations of the susceptibihty 
and specific- heat peak maximum to the infinite volume limit yielded [0], 
respectively, Pc = 0.6930(2) and (3c = 0.6931(10), in agreement with the 
cumulant crossing value. 

However, PFL did not mention a later reanalysis of the high-tem- 
perature series expansion based on the Pade {Pc = 0.6924(2)) and ratio 
{(3c = 0.6925(1)) method, respectively. While these values are consistent 
with the TM estimates given in ref.0, we are now faced with the problem 
that the latest MC result is significantly higher. 

The critical coupling is a non-universal parameter and from this point of 
view not of particular interest. Most estimates of universal critical indices, 
however, are biased and usually depend quite strongly on the precise value 
of (3c- To clarify the above discrepancy we found it therefore worthwhile to 
perform an independent high precision MC study on larger lattices of sizes 
up to 48^ and with even higher statistics as in PFL's work. Relying on the 
Metropolis algorithm as PFL did, such a project would have been hardly 
feasible. As with most local algorithms, the Metropolis (pseudo) dynamics 
suffers from the severe problem of critical slowing down, that is from large 
autocorrelation times r = aL^ (with dynamical critical exponent z = 1.94(6) 
(ref.0) and[| a ^ 3.76), which reduce the size of the statistical sample, N, 
effectively to N^.^ = N/2t. It is therefore crucial to use one of the improved 
algorithms of the past few years that avoid this problem. We chose the 

^ This can be read off from Fig. 2 in ref. M . 
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cluster algorithm |]rT| in its single-cluster variant |T^. From studies of related 
spin models it is known [1^ that this update algorithm is extremely efficient 
in three dimensions. 

The classical Heisenberg model is defined by the partition function 



n 



id cos 6i 



where /3 = 1/T is the (reduced) inverse temperature and 



(1) 



(2) 



is the total energy (in the following small letters always denote the corre- 
sponding intensive quantities, e.g., e = E/V). The sum runs over all nearest 
neighbour pairs and the three-dimensional unit spins s at the sites i of 
a simple cubic lattice are parametrized as s = (cos sin 6', sin sin 6', cos 6*). 
We always employ periodic boundary conditions. 

Our simulations were organized as follows. First, we did one run for 
each lattice size at /5o = 0.6929, the estimate of /3c by PFL, and recorded 
the energy histogram Pj3g{E) and the microcanonical averages {{m'')){E) = 
J2m Ppq{E, M)m^ / Pp^{E) , k = 1,2, 4, where m = \rn\ is the magnitude of the 
magnetization m = ^ J2x ^(^) of ^ single spin configuration. The tempera- 
ture independent averages {{m^)){E) can be computed by accumulating the 
values of in lists indexed by the associated energy bin of the configuration 
and normalizing at the end by the total number of entries in each bin, making 
it thus unnecessary to store the two-dimensional histogram Pj3g{E,M). The 
continuous energy range < E < 3V was discretized into 90000 bins. The 
data of this run is sufficient to compute the approximate positions /3_ < j3o 
and /5+ > /?o of the (connected) susceptibility and the specific-heat peak 
maximum by reweighting techniques . We then performed two more runs 
at /5_ and respectively, again recording Ppi^E) and {{m'')){E). This choice 
has the advantage that one automatically stays in the critical region since 
both P_ and /3+ scale with L~^/^, where u is the correlation length exponent. 

-,0, + for all 



From this data we can compute three estimates 0^i\j3) 
thermodynamic observables Ol of interest, and for any (3 value in the vicin- 
ity of Po, /3+ by reweighting . Furthermore, since we devided the whole 

run into several blocks and stored the energy histograms and microcanon- 

)(") 
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ical averages for each block, we could compute jackknife errors AO)^ 
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on 0^^\ This allowed us to get an optimized average of these three values 
that minimizes the relative error of the combined Ol{(3) for each observable 
separately (the relative weights are simply l/(AO["^)2). All our runs contain 
at least 10000 x r measurements, where r is the integrated autocorrelation 
time of the susceptibility. As expected for the single-cluster update, r turns 
out to be almost independent of the lattice size and to be very small (< 2, 
in units of lattice sweeps that allow direct comparison with the Metropolis 
algorithm). For the 48'^ lattice, r is about three orders of magnitude smaller 
than for the Metropolis algorithm. This explains why we could study much 
larger lattice sizes than PFL, and could still afford to have about ten times 
better statistics. 

To determine (3c we first concentrated on Binder's cumulant 

t'.W-l-^. (3) 

where the angular brackets (. . .) denote thermal averages with respect to ([l|). 
It is well known that, asymptotically for large L, all curves Ul{(3) should 
cross in the unique point {[3c, U*). Our results for lattices of size L=12, 16, 
20, 24, 32, 40, and 48, obtained by the optimization procedure described 
above, are shown in Fig. 1. We see that all curves indeed cross each other at 
approximately the same (3 value of 0.693. A closer look, however, reveals that 
the crossing points of f/z. and U l' are systematically shifted, depending on the 
ratio b = L/L'. This is the expected behaviour for finite lattices, caused by 
confluent correction terms. Employing well-known [0 extrapolation formulas 
we obtain the final estimates |]15| 

Pc = 0.6930 ± 0.0001, (4) 

and 

[/* = 0.6217 ±0.0008. (5) 

The critical coupling is thus found in excellent agreement with the value 
quoted by PFL, and also U* agrees very well with their estimate of 0.622(1). 
For comparison, a field theoretic expansion in y/A — D predicts for D = 3 a 
4% lower value of 0.59684... 



Let us now turn to the FSS estimates of critical exponents. The deriva- 
tives dUL/d(3 at Pc = 0.6930 should scale asymptotically for large L with 
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L^/'^. In a log- log plot of all our data points we find a perfect straight line 
fit (with goodness-of-fit parameter Q = 0.61) and from the inverse slope we 
read off 

z/ = 0.704 ± 0.006, (6) 

which again is in agreement with the value quoted by PFL, u = 0.706(9) 
(determined by the same method, but at /? = 0.6929), and with the field 
theoretical estimates of z/ = 0.705(3) (resummed perturbation expansion 
[01 )) ^ = 0.710(7) (resummed e-expansion ||l^). The high quality of this 



fit (as well as of all other fits described below) shows that the asymptotic 
scaling formula works down to our smallest lattice size L = 12, indicating 
that there is no need for confiuent correction terms. 

The ratio of exponents f3/iy follows from the scaling of the magnetization, 
(m) oc L^^^'^. In a log-log plot of (m) at Pc = 0.6930 vs L we obtain from a 
straight line fit (with Q = 0.68) the estimate 

= 0.514 ± 0.001, (7) 

which is slightly lower than the value given by PFL, jS/u = 0.516(3) (deter- 
mined at /5 = 0.6929). To test by how much our result is biased by the value 
of (3c we have redone our analysis at /? = 0.6929. Here we get the slightly 
higher value of 0.519(1). The quality of the fit, however, is much worse 
{Q = 0.30). Since we observe a similar worsening of the fit at /? = 0.6931 
(/3/i/ = 0.509(1), Q = 0.31), we take this as support for our estimate of Pc- 
We rely on the goodness-of-fit parameter since visually it is impossible to 
make a distinction between these fits when plotted on a natural scale. It 
should be emphasized that even these slight variations in the estimate of the 
critical coupling produce significant changes in the estimate of the exponent 
ratio that clearly dominate the statistical errors. Combining @ and (0) we 
get for the critical exponent (3 = 0.362(4), to be compared with the estimates 



of /5 = 0.3645(25) (ref.^) and f3 = 0.368(4) {leL^ 

All other critical exponents follow from (hyper-) scaling relations. In this 
way we obtain 

a/z/ = 2/z/ - D = -0.159 ± 0.025, (8) 
7/z/ = 3-2/3/z/= 1.972 ±0.002, (9) 

and 

= 2 - 7/z/ = 2/3/z/ - 1 = 0.028 ± 0.002. (10) 
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To test these relations we have performed also a direct analysis of the 
critical behaviour of the specific heat 

C = V(3^ ((e^) - (e)2) oc L"/^ (11) 

and of the (connected) susceptibility 

= V(3 ((m^) - (m)2) oc L^^". (12) 

Since the magnetization vanishes at and above the critical temperature we 
may also take 

X = V(3{m^) oc L^/^ (13) 

which is usually a less noisy estimator. 

From the location of the maxima of C and we can get further estimates 
for the critical coupling by assuming the FSS relation Tmax = Tc + aL~^/'^ + .... 
Using our value of u = 0.704 we obtain from the linear fits shown in Fig. 2 
the estimates = 0.6925(9) (from Tc^^^ with Q = 0.80) and = 0.6930(3) 
(from T^5^^^ with Q = 1.0), respectively. These values are consistent with the 
cumulant crossing value (^, but have larger statistical errors. 

From the scaling of x'^ and x at our best estimate of (3c = 0.6930 we get 
from linear fits t] = 0.0156(44) (Q = 0.69) and r] = 0.0271(17) (Q = 0.78), 
respectively. The latter value is in perfect agreement with the scaling pre- 
diction ([10|). Moreover, performing fits to x at /? = 0.6929 and /? = 0.6931 
we obtain estimates of = 0.0364(17) (Q = 0.36) and r] = 0.01781(17) 
{Q = 0.43), respectively, that are again in perfect agreement with the scaling 
predictions based on the corresponding fits to the magnetization. This is 
not unexpected since the measurements of (m) and x are of course strongly 
correlated. Finally, analyzing the FSS behaviour of the susceptibility maxi- 
mum, Xmax L'^^" , we estimate rf = 0.0231(61) {Q = 0.60). Notice that all 
MC estimates are lower than the field theory values which are rj = 0.033(4) 
(ref. [0) and r] = 0.040(3) (ref. [|T§)g. 

The specific heat exhibits a finite, cusp-like singularity, because a has 
a negative value. We tried a three-parameter fit of the form Cmax = a — 

^ In a recent recalculation ||l9| of all Feynman graphs several errors in the highest order 
of the e-expansion were corrected. A subsequent reanalysis of the resummed series 
gave a slightly smaller value for rj. Compared with the error bar, however, this correction 
is negligible. 
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jjl^a/u^ The result a/u = —0.33(22) (Q = 0.69) is compatible with the scaling 
relation but, due to its very large statistical error, it does not provide a 
stringent test of this scaling prediction. Another way of testing eq. is to 
assume the predicted value of a/u = —0.159 and to fit only the parameters 
a and b. The resulting fit turned out to be of almost equally good quality 

In summary, using high-precision data from single-cluster MC simulations 
combined with optimized multi-histogram techniques, we have performed a 
fairly detailed FSS analysis of the 3D classical Heisenberg model on simple 
cubic lattices of size up to 48^. Qualitatively, our main result is that the 
asymptotic FSS region sets in for small lattices sizes, L ?a 12. Quantita- 
tively, our value for the critical couphng, /3c = 0.6930(1), is in almost perfect 
agreement with the MC estimate reported recently by PEL (ref. 0), but 
is significantly higher than estimates from high-temperature series expan- 
sion analyses and transfer matrix methods. Our results for the two basic 
critical exponents, u = 0.704(6) and j3 = 0.362(4), are in good agreement 
with field theoretic predictions. Scaling relations imply a = —0.112(18) and 
rj = 2 — 'J /u = 0.028(2). Direct measurements of these exponents provide 
tests of the scaling relations. In the case of rj we find good agreement when 
the scaling of x a-t Pc is considered. Using at Pc or Xmax^ however, the 
situation is less clear. In the case of a, its negative value causes numeri- 
cal problems, since a finite, cusp-like singularity is notoriously difficult to 
analyze. 
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Figure Headings 



Fig. 1: Fourth-order cumulant Ul vs p. The values of Ul{P) were obtained 
by reweighting and optimized combining the results of our three simula- 
tions at different temperatures for each lattice size L. The simulations 
were performed at Po = 0.6929 (the critical coupling found by Peczak et 
al. [^), and at the temperature locations of the maxima of the specific 
heat C and the susceptibility respectively. 

Fig. 2: Variation of the pseudo transition temperatures T^f^^ (L) and Tc,„^^ (L) 
with L~^/^, where u = 0.704(6) is our FSS estimate (see text). The 
fits yield estimates of (3^ = 0.6930(3) (Q = 1.0) and f3c = 0.6925(9) 
{Q = 0.80), respectively. 
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